% Barcode Recognition example 
% ----------------------------
%
% This is the recognizer for Interleaved 2 of 5 barcode. 
%
% The recognizer reads the image from file system, fixed location. The image is expected
% to contain correctly rotated image of a barcode, the horizontal centre of the image must
% contain the barcode.  No false bars are expected to be present, the algorithm does not
% detect fals barcodes, it assumes the detected bars belong exclusively to barcode.
%
% The recognition strategy is optimized for mobile devices. I uses a simple threasholding
% strategy to deliver binary image, detects the bars, determines the structure and decodes
% the aplhabet.
%
% 	1. Image is loaded into memory. Algorithm finds the horizontal mid section, extract
% 	   single row of signal and determines the threashold by using three different methods:
%
% 	   a. calculates the mean value of the signal (all colors) and determines.
%
% 	   b. finds all the minimums and maximum in the signals (peakdet), determines the
% 	      avegare for both maximums and minimums, and aritmetics average for averaged
% 	      maximums and averaged minimums.
%
% 	   c. Finds the maximum in the signal, finds the minimum in the signal and determines
% 	      arithmetics average for it.
%
% 	   All the mentioned three methods are used to binarize the image (single row of image
%		of image data) and the binarized image if further processed as follows:  
%
% 	2.	The bars are detected. The bars detection uses simple XOR function and records bar
% 	   every time the binary signal chnages polarity.
%
% 	3. 
%
% (used test image from iPhone worked fine for broad interval of binarization thresholds
		% ranging from 0.25 - 0.42 (64-108) while the mid level threashold (128) delivers
		% very noisy data)
%
% This example performs very simple binarization and single line extraction. The barcode
% is then detected on single line only.  Bars are extracted from the filterred image. The
% leading white space is removed, barcode leading sequecne detected and removed.
% Frequency analysis is applied to the signal fo the bars. The first maximum corresponds
% with the width of the narrow bar. The frequency data is fitted by a sine function and
% the model is searched for local maximum. Its position yields the ideal width of the
% narrow bar.  Barcode is then normalized and weighted, symbology determined and symbols
% converted into decimal digits based on the barcode symbology.



clear;

% Aplhabet of Interleaved 2 of 5
alphabet = [
	1,1,2,2,1;
	2,1,1,1,2;
	1,2,1,1,2;
	2,2,1,1,1;
	1,1,2,1,2;
	2,1,2,1,1;
	1,2,2,1,1;
	1,1,1,2,2;
	2,1,1,2,1;
	1,2,1,2,1];

%<<<
function err = drop_bc()
	printf("drop_bc\n");
	clear X;
	clear lf;
	clear res;
	clear bars;
	clear nbars;
	clear fin;
	clear lfin;
	clear nlfin;
	clear y;
	clear x;
	clear x1;
	clear x2;
	clear a;
	clear yc;
	clear ys;
	clear bc;
	clear lead;
endfunction
%>>>

figures = input("Show figures? [0|1,0|1]: " );

% Read image from FS ...
%I=imread("./code128_ios.jpg");
I=imread("./barcode_0.jpg");
disp("Image size:");
disp(size(I));

%Identify the cut through the barcode ...
cut=uint16(rows(I)/2);

detected = false;
while( not(detected) )
	cut = cut + 8;
	disp("Selected cut through the barcode:");
	disp(cut);
	if( uint16(cut) > uint16(rows(I)) )
		break;
	endif

	% Extract color information from the image
	sr=I(cut,:,1);
	sg=I(cut,:,2);
	sb=I(cut,:,3);
	
	% Determine the binarization factor by finding the mid 
	% point for the intensity signal for all three colors.
	th0 = uint16(max(sr))+uint16(max(sg))+uint16(max(sb));
	ths(1,3) = single(th0)/(6*256);
	
	disp("Suggested threashold intensity mid point of max:");
	disp(ths(1,3));
	
	th0 = uint16(mean(sr)) + uint16(mean(sg)) + uint16(mean(sb));
	ths(1,1) = single(th0)/(3*256);
	
	disp("Suggested threashold intensity based on mean: ");
	disp(ths(1,1));
	
	[maxs,mins]=peakdet(int16(sr),60);
	top=mean(maxs(:,2));
	bottom=mean(mins(:,2));
	ths(1,2)=(top-bottom)/512;

	disp("Suggested threashold intensity based on peakdet: ");
	disp(ths(1,2));

	% Create three binarized images based on three threasholds ... 
	for thi=1:1
%		th = th1 - (0.25*(th1-th0));
		th = ths(1,thi);
		disp("Binarizing the image with threashold intensity : ");
		disp(th);
		
		%binarize the image with the pre-set threashold.
		BW=im2bw(I,th);
		disp("Binarized image size:");
		disp(size(BW));
		
		%identify the position of the mid cut through the barcode
		l=BW(cut,:);
		
		disp("size of barcode area:")
		disp(size(l));
		
		%save barcode strip for optical check
		%imwrite(l,"/tmp/exp.tiff","tif");
		%disp("Image saved into /tmp/exp.tiff");
		
		%despecle the image by removing 1x1 speckels, both white and black
		for i=2:columns(l)-1
			lf(1,i) = despecle(l(i-1),l(i),l(i+1));
		endfor
		
		%pad the image
		lf(1,1) = lf(1,2);
		lf(1,columns(l)) = lf(1,columns(l)-1);
		
		%save the filter result for optical check
		imwrite(lf,"/tmp/expf.tiff","tif");
		disp("Image saved into /tmp/expf.tiff");
		
		% swap the filtered and original image
		bwi(thi,:) = lf;
	endfor

	l = bwi(1,:);
	
	% Draw the binarization step result on the figure ....
	if( figures(1,2) ) %<<<
		figure;
		subplot(2,1,1);
		plot(uint16(1:columns(sr)),uint16(sr(1,1:end)),"-r");
		hold on;
		plot(uint16(1:columns(sg)),uint16(sg(1,1:end)),"-g");
		hold on;
		plot(uint16(1:columns(sb)),uint16(sb(1,1:end)),"-b");
		xlabel("Image X-coordinate [px]");
		ylabel("Image RBG signal");
		legend("analog RGB signal");
		subplot(2,1,2);
		bar(uint16(l));
		xlabel("Image X-coordinate [px]");
		ylabel("Image binary signal");
		legend("binarized inverted signal");
	endif
   %>>>
	
	%Identify bars by running edge detection algorithm. Build
	% the signal for bars as a response signal for edge detection
	pos = 1;
	p = 1;
	for i=2:columns(l)
		if( bitxor(l(i-1),l(i)))
		 	res(pos) = p;
			p = 1;
			pos = pos + 1;
		else
			p = p + 1;
		endif
	endfor
	
	%strip first bar as the leading blank space.
	% ideally the width is checked, not in this case :-(
	bars=res(1,2:end);
	
	%normalize for minimal width
	disp("Minimum width=");
	disp(min(bars));
	disp("Maximum width=");
	disp(max(bars));
	
	% Calculate the frequency spectrum for detected bars
	fin(1,1:max(bars)) = 0;
	for i=1:columns(bars)
		fin(bars(1,i)) = fin(bars(1,i)) + 1;
	endfor
	
	% Draw the bar detection results ...
	if( figures(1,1) )
		figure;
		subplot(2,1,1);
		bar(bars);
		ylabel("Amplitude [px.]");
		xlabel("Bar Count [No.]");
		title("Barcode signal and its spectrum");
		legend("bars");
		subplot(2,1,2);
		bar(fin);
		ylabel("Frequency [No.]");
		xlabel("Bar Width [px]");
		legend("bars");
	endif
	

	% the model is masked against the trimmed mask. Trimming of the mask 
	% uses normalized vector and removes every bar with frequency bellow 10%.
	% First after the trimming and filtering the model is searched for maximum.
	lfin=fin(1,1:(columns(fin)/2));
	nlfin = lfin/(max(lfin));
	for i=1:columns(lfin)
	   if( nlfin(1,i) > 0.1 )
			nlfin(1,i) = 1;
		else
			nlfin(1,i) = 0;
		endif
	endfor
	

	% Once the mask is trimmed, the filter runs and removes specles (1x1) in mask.
	% Specle means the missing bar (single bar with frequency less than 10%
	% is not removed, is considered present).
	for i=2:columns(nlfin)-1
		nlfin(1,i) = despecle(nlfin(i-1),nlfin(i),nlfin(i+1));
	endfor
	nlfin(1,1) = nlfin(1,2);
	nlfin(1,columns(nlfin)) = nlfin(1,columns(nlfin)-1);
	
	ye = regsin0( fin );
	yc = regsin1( fin );	
	ys = regsq( fin );	

	% finally the mask is run against the model which removes all the side lobes
	% produced by fitting model. It is done by eleme by element masking.
	nyc=yc.*(nlfin'); % sin(x) model
	nys=ys.*(nlfin'); % quadratic model
	nye=ye.*(nlfin'); % sin(x) model
	
	% This is important.
	% Sin(x) can be phase shifted (pi), and the interpolation actually
	% shows minimum where the maximum shall occur. This algoriths reveales
	% the phase shift and swaps the minimums and maximums.
	[maxs,mins]=peakdet(nyc,0.5);
	if( rows(maxs) > rows(mins) && rows(maxs) < 3 )
		disp("Phase shifted sin(x) encountered");
		narrow_bar_widths(1,1) = mins(1,1);
	else 
		disp("Phase correct sin(x) encountered");
      if( rows(maxs) == 3 )
   		narrow_bar_widths(1,1) = maxs(2,1);
      else
         for j=1:maxs(:,1)
            if( maxs(j,2) == max(maxs(:,2)))
         		narrow_bar_widths(1,1) = maxs(j,i);
            endif
         endfor
      endif
	end
	
	% the masked model data is searched for position of the maximum,
	% The position corresponds width the ideal width of the narrow bar.
	% for sin model
	%{
	for i=1:rows(nyc)
		if( nyc(i) == max(nyc))
			narrow_bar_widths(1,1) = i;
		end
	end
	%}

	for i=1:rows(nye)
		if( nye(i) == max(nye))
			narrow_bar_widths(1,2) = i;
		end
	end
	
	% ... and for quadratic model
	for i=1:rows(nys)
		if( nys(i) == max(nys))
			narrow_bar_widths(1,3) = i;
		end
	end
	
	% Do the representative drawings ...
	disp("Detected width of the narrow bar improved sin(x)");
	disp(narrow_bar_widths(1,1));
	disp("Detected width of the narrow bar old sin(x)");
	disp(narrow_bar_widths(1,2));
	disp("Detected width of the narrow bar x^2=");
	disp(narrow_bar_widths(1,3));
	
	if( figures(1,1) )
		figure;
		subplot(1,2,1);
		bar(lfin,"facecolor", "g", "edgecolor", "b");
		hold on;
		plot(yc,"-b");
		plot(ys,"-m");
		plot(ye,"-r");
		ylabel("Frequency [No.]");
		xlabel("Bar Width [px]");
		title("Frequency profile for narrow bar");
		legend("data", "improved sin(x)", "polynom", "sin(x)");
		subplot(1,2,2);
		bar(nlfin,"facecolor", "g", "edgecolor", "b");
		hold on;
		plot(nyc,"-b");
		plot(nys,"-m");
		plot(nye,"-r");
		ylabel("Frequency [No.]");
		xlabel("Bar Width [px]");
		title("Normalised Frequency profile for narrow bar");
		legend("mask", "improved sin(x)", "polynom", "sin(x)");
	endif
	
	% First try toi quadratic result agaist the given binarization
	% Assuming the barcode is not detected, then change the width
	% to data recommended by the sin(x) model. Bails out if
	% none of the models deliver complete barcode.
	
	narrow_bar_width=narrow_bar_widths(1,1);
	narrow_bar_width_index = 1;
	
	disp( "Trying to identify barcode, using the bar width:");
	disp( narrow_bar_width );
	
	do 
		% Vector of the detected bars is normalized against the
		% ideal width and rounded to nearest integer.
		nbars=bars/narrow_bar_width;
		nbars=round(nbars);
		
		%separate the barcode signature from digits
		bc=nbars(1,5:end);
		lead=nbars(1,1:4);
		disp("Leading sequence:");
		disp(lead);
		
		if( lead == [1,1,1,1] ) 
		   disp( "Barcode Interleaved 2of5 detected");
			% The barcode vector is normalized for narrow (1) and wide (2) bars. Thsi is OK for 
	      % all barcodes that uses just two widths of the bars (narrow and thick).
			for i=1:columns(bc)
				if( bc(i) < 1 )
					bc(i) = 1;
				endif
				if( bc(i) > 1 )
					bc(i) = 2;
				endif
			endfor
			
			%decode the bars into barcode symbols. Save barcode symbols in 'd' matrix
			pos = 1;
			for i=1:10:columns(bc)-10
				v=[bc(1,i),bc(1,i+2),bc(1,i+4),bc(1,i+6),bc(1,i+8)];
				w=[bc(1,i+1),bc(1,i+3),bc(1,i+5),bc(1,i+7),bc(1,i+9)];
				d(pos,:)=v;
				d(pos+1,:)=w;
				pos = pos + 2;
			endfor
			
			%decode the barcode symbols into decadic numbers
			for i=1:rows(d)
				BC(i) = -1;
				for j=1:10
					if( d(i,:) == alphabet(j,:) ) % Vector aritmetics
						BC(i) = (j-1);
					endif
				endfor
			endfor
	
	      detected = true;

		   %
		   % Show results!
		   disp(BC);
	
	      % Check whether all the digits got recognized, if not use the other digit.
	      if( min(BC) == -1 )
	         detected = false;
				narrow_bar_width_index = narrow_bar_width_index + 1;
	         if( narrow_bar_width_index == 4 )
		      	disp( "next cut" );
					drop_bc();
					narrow_bar_width_index = 1;
	         	break;
				endif;
	         narrow_bar_width=narrow_bar_widths(1,narrow_bar_width_index);
	         disp("Barcode not recognized successfully. Changing bar width");
	         disp( narrow_bar_width );
	      endif
	
		else
	      detected = false;
		   disp( "Barcode Interleaved 2of5 not detected at all");
			narrow_bar_width_index = narrow_bar_width_index + 1;
	      if( narrow_bar_width_index == 4 )
		      disp( "Next cut ..." );
				drop_bc();
				narrow_bar_width_index = 1;
	         break;
			endif;
	      narrow_bar_width=narrow_bar_widths(1,narrow_bar_width_index);
		   disp( "Changing the bar width to:");
	      disp( narrow_bar_width );
		endif
	until( detected );

endwhile
disp("*STOP*");

